% This file generates Figure B1 in Appendix B of the paper. It presents the
% IRF as a function of the duration of ZLB from a Taylor Rule.

close all; clc; clear;

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Parameters and Steady-State Values %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


beta                    =   0.99;   % Discount factor
chi                     =   0;      % Inverse of Frisch elasticity. 
sigma                   =   1;      % Inverse of EIS.

phi                     =   0.75;   % Prob. firms cannot adjust prices. 

rho                     =   0.9;    % Persistence of natural rate shock. 

% Slope Phillips Curve:
gamma                   =   ((1-phi)*(1-phi*beta)/phi)*(sigma+chi);

alpha                   =   1.5;    % Coefficient of inflation in the rule. 
% Auxiliar equation:
delta                   =   sigma*(1+chi)*(rho-1)/(sigma+chi);

% Output gap policy:
theta_1                 =   ((1-beta*rho)/sigma)/((1-rho)*(1-beta*rho)...
                            +gamma/sigma*(alpha-rho))*delta;
% Output gap policy:
theta_2                 =   (gamma/sigma)/((1-rho)*(1-beta*rho)...
                            +gamma/sigma*(alpha-rho))*delta;

T                       =   21+1;   % Number of periods.           

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             Creating IRFs                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

H                       =   3;      % Peg length. 

% Pre-allocating:
pi                      =   ones(H+1,1);
x                       =   ones(H+1,1);
TFP                     =   ones(T,1);
infl_1                  =   ones(T,1);
output_1                =   ones(T,1);
% Initial values:
x(H+1)                  =   theta_1*rho^H;
pi(H+1)                 =   theta_2*rho^H;
a                       =   1;

% Creating IRFs:
for ix = H:-1:1                     % We go backwards. 
    x(ix)                   =   x(ix+1)+1/sigma*(pi(ix+1)+delta*rho^(ix-1));
    pi(ix)                  =   gamma*x(ix)+beta*pi(ix+1);
end

for ix = 1:T
    TFP(ix)                 =   rho^(ix-1)*a;
    if ix < H+1
        infl_1(ix)              =   pi(ix);
        output_1(ix)            =   x(ix)+(1+chi)/(sigma+chi)*TFP(ix);
    else
        infl_1(ix)              =   theta_2*TFP(ix);
        output_1(ix)            =   theta_1*TFP(ix)+(1+chi)/(sigma+chi)*TFP(ix);
    end
end

%% Different Peg. 
H                       =   6;

% Pre-allocating:
pi                      =   ones(H+1,1);
x                       =   ones(H+1,1);
TFP                     =   ones(T,1);
infl_2                  =   ones(T,1);
output_2                =   ones(T,1);
% Initial values:
x(H+1)                  =   theta_1*rho^H;
pi(H+1)                 =   theta_2*rho^H;

% Creating IRFs:
for ix = H:-1:1
    x(ix)                   =   x(ix+1)+1/sigma*(pi(ix+1)+delta*rho^(ix-1));
    pi(ix)                  =   gamma*x(ix)+beta*pi(ix+1);
end

for ix = 1:T
    TFP(ix)                 =   rho^(ix-1)*a;
    if ix < H+1
        infl_2(ix)              =   pi(ix);
        output_2(ix)            =   x(ix)+(1+chi)/(sigma+chi)*TFP(ix);
    else
        infl_2(ix)              =   theta_2*TFP(ix);
        output_2(ix)            =   theta_1*TFP(ix)+(1+chi)/(sigma+chi)*TFP(ix);
    end
end

%% Different Peg. 
H                       =   10;

% Pre-allocating:
pi                      =   ones(H+1,1);
x                       =   ones(H+1,1);
TFP                     =   ones(T,1);
infl_3                  =   ones(T,1);
output_3                =   ones(T,1);
infl_flex               =   ones(T,1);
output_flex             =   ones(T,1);
% Initial values:
a                       =   1;
x(H+1)                  =   theta_1*rho^H;
pi(H+1)                 =   theta_2*rho^H;

% Creating IRFs:
for ix = H:-1:1
    x(ix)                   =   x(ix+1)+1/sigma*(pi(ix+1)+delta*rho^(ix-1));
    pi(ix)                  =   gamma*x(ix)+beta*pi(ix+1);
end

for ix = 1:T
    TFP(ix)             =   rho^(ix-1)*a;
    if ix < H+1
        infl_3(ix)          =   pi(ix);
        output_3(ix)        =   x(ix)+(1+chi)/(sigma+chi)*TFP(ix);
    else
        infl_3(ix)          =   theta_2*TFP(ix);
        output_3(ix)        =   theta_1*TFP(ix)+(1+chi)/(sigma+chi)*TFP(ix);
    end
end

for ix = 1:T
    infl_flex(ix)       =   theta_2*TFP(ix);
	output_flex(ix)     =   theta_1*TFP(ix)+(1+chi)/(sigma+chi)*TFP(ix);
end

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           Plotting results                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

t                       =   0:T-2;

figure
subplot(1,2,1)
plot(t,output_flex(1:T-1),'-b',t,output_1(1:T-1),'--b',t,output_2(1:T-1),':b',t,output_3(1:T-1),'-.b','Linewidth',1.5)
xlabel('Horizon','Interpreter','latex','fontSize',14)
ylabel('Output','Interpreter','latex','fontSize',14)
h = legend('H=0 (No ZLB)','H=3','H=6','H=10','Location','SouthEast');
set(h,'Interpreter','latex')
grid on

subplot(1,2,2)
plot(t,infl_flex(2:T),'-b',t,infl_1(2:T),'--b',t,infl_2(2:T),':b',t,infl_3(2:T),'-.b','Linewidth',1.5)
xlabel('Horizon','Interpreter','latex','fontSize',14)
ylabel('Expected Inflation','Interpreter','latex','fontSize',14)
h = legend('H=0 (No ZLB)','H=3','H=6','H=10','Location','SouthEast');
set(h,'Interpreter','latex')
grid on
set(gca,'layer','top')
hold

